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ABSTRACT 


Possibly the best set of data for earthquake excitation of soils exists for the test site operated 
by the Taiwan Power Company in conjunction with the Electric Power Research Institute (EPRI) 
at Lotung Taiwan. At this site, two locations are instrumented with three-component 
accelerometers at depths of 47 m, 17 m, 11 m, 6 m, and at the surface. One array is in the free- 
field, while the other is adjacent to a one-quarter scale nuclear containment vessel. The site is also 
well instrumented with piezometers at various depths and locations. The simplified soil profile 
consists of 30 m to 35 m of silty sand and sandy silt with some gravel, overlaying a thick clay and 
silt deposit. The water table is within half a meter of the ground surface. This area is seismically 
active, and strong shaking generated by many earthquakes exhibiting a wide range of magnitudes 


have been recorded since 1986. 


For this study the modal frequencies and damping ratios were calculated for events 3, 4, 7, 
8, 9, 10, 12 and 16 with local magnitudes ranging from 4.5 to 7.0. The modal frequencies and 
damping ratios calculated are examined for the effect of local energy intensity and soil structure 
interaction. Modal frequencies are seen to decrease with increasing intensity once a certain 
threshold of acceleration/intensity is reached. This result is consistent with the data obtained by 
other authors using different techniques. For the 0-6 m interval the decrease in frequency with 
event energy is less pronounced under a model containment structure than in the free field. This 
soil-structure effect is increasingly diminished with depth and absent by the 17-47m interval. 
Calculated damping values demonstrate an expected increase with input seismic energy. For the 
0-6 m and 6-11 m intervals the damping values are higher under the model structure than in the 
free field. This distinction is completely missing in the 17-47 m results. The transition to non- 
linear behavior, while less pronounced with increasing depth, consistently occurs above a peak 


acceleration of 0.05 g or Anas Intensity of 100 m/sec. 


The results clearly indicate a degree of non-linear response over the intervals studied. 
Evidence of a decrease in specific interval fundamental frequency and an accompanying general 
trend of increased damping with higher seismic energy are clear. Comparison of the results of this 
study with previous work considered with the inherent superiority of parametric modeling for 
transient and/or non-stationary time series such as earthquakes indicate that system identification 
is a more robust method for identifying fundamental frequencies and damping values for layers of 


earth materials when borehole information is available. 
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CHAPTER 1 - INTRODUCTION 


1.1 Background 


There are many unanswered questions of interest to the geotechnical community 
concerning the behavior of soils subjected to earthquake excitation. Foremost among them are 
questions concerning the strain-dependent non-linear behavior of soils, and soil-structure 
interaction. In an attempt to gain further rational understanding of these problems the Electric 
Power Research Institute (EPRI), the Colorado School of Mines (CSM), and the Structures 
Division of the National Institute of Standards and Technology (NIST) formed a cooperative 


research team to evaluate ground motion time histories recorded at the Lotung site, Taiwan. 


Much of the necessary data was made available by EPRI from the Lotung site. The Lotung 
strong motion data set are an extremely unique set of data. The completeness of this input-output 
data set makes it ideal for analysis using system identification (SI) methods. Data from the 
Wildlife Site, Imperial Valley CA, and Treasure Island are freely available. This report lays out 
the work undertaken with funding provided by NIST under Award Number 60NANB500074. 


1.2 Why Use System Identification? 


An important goal for earthquake engineering is the ability to estimate soil properties 
without intruding into the soil mass. For the engineer interested in seismic behavior of soils, the 
dynamic properties of the soil are of interest, particularly large-strain properties. The archetypal 
large strain field excitation is earthquake strong motion. Ideally, both ground motions into the soil 
layer of interest and on the surface above the layer are recorded, as illustrated by the cartoon in Fig. 
1.1. Given this known input propagating upward from depth, and the output at the top of the soil 
column, the behavior of the soil can be modeled by inverse theory. If a suitable model is chosen 
to represent the system of interest, the estimated model parameters will correspond to important 
mechanical parameters of the system, such as damping, natural frequency, and stiffness. This 


estimation of parameters is commonly known as system identification (SI). 


Output receiver 


WORN GIT. 
System Parameters to be Determined 
damping 
frequency 
stiffness 


de receiver 
Earthquake shaking 
Fig. 1.1 - Configuration of the System Identification Method. 


The traditional method of geotechnical analysis of dynamic soil motions is through the 
Fourier transform. However, serious problems arise when this method is applied to short data 
streams, and to signals changing through time — non-stationary signals. This study was 
undertaken to show the effectiveness of a different type of model, a parametric model commonly 
used in automatic control and geophysics, which avoids many of the limitations inherent with 
calculating the system transfer function by Fourier techniques. An important aspect of this 
particular parametric model is the theoretical link between the estimated system parameters and 
the mechanical parameters of a lumped-mass oscillator. The parametric model allows estimates of 


system dynamic properties to be made if an input-output data set is available. 
1.3 Scope 


This report presents detailed SI analyses of the full suite of seismic events at the Lotung 
site. Estimates of resonant frequencies and damping, and the dependence of these quantities on 
earthquake intensity will be presented. To this end the Lotung site itself will be introduced in 
Chapter 2 through geological, seismological, and geotechnical descriptions. The SI procedures 
will be discussed in detail in Chapter 3. Chapter 4 presents results of the system identification of 
the Lotung site data, with a discussion and conclusions given in Chapter 5. Appendix A describes 
the signal processing used on the raw data, with the integration procedure given in Appendix B. 
The pore water pressure time histories are shown in Appendix C. A literature review of damping 
values measured by other researchers is compiled in Appendix D. Finally, comparisons of actual 


interval outputs and modeled outputs for all events and intervals is given in Appendix E. 


CHAPTER 2 - THE LOTUNG SITE, TAIWAN 


pats | Introduction and Geography 


With the growth of the use of nuclear-powered generating plants in the 1970’s, many safety 
related questions about the seismic performance of these plants arose. In the early 1980’s, EPRI 
and the Taiwan Power Co. constructed two scale models (1/4 and 1/12 scale) of a nuclear 
containment structure near Lotung, Taiwan. This is a very seismically active area in northeast 
Taiwan (see Fig. 2.1). The site and structures were elaborately instrumented so that soil and 
structural response, and soil-structure interaction, to earthquakes could be carefully studied (Tang 


et al., 1989; Liu and Yeh, 1985). 
pe Instrumentation 


The soil instrumentation includes a three-arm surface array, as shown in Fig. 2.2a. The 
arms radiate approximately 47 m from the 1/4 scale containment structure. In addition, there are 
two downhole arrays of accelerometers extending to a depth of 47 m, as shown in Fig. 2.2b. The 
surface accelerometers are triaxial force-balance units (Kinemetrics FBA-13) oriented in the N-S, 
E-W, and vertical directions. The downhole arrays (DHA and DHB) are modified Kinemetrics 
FBA-13H units oriented in the N-S, E-W, and vertical directions. DHA is located 3 m from the 
containment vessel and DHB is located 47 m from the structure, allowing identification of the 
effects of the structure on soil response. The downhole instruments are located at depths of 6 m, 
11 m, 17 m, and 47 m. The simplified soil profile consists of 30-35 m of silty sand and sandy silt 
with some gravel, above clayey silt and silty clay. The water table is within half a meter of the 


ground surface, under artesian pressure. 
2.3 Site Characterization 


The basic geology of the Lotung site is summarized by Wen and Yeh (1984) and Tang 
(1987). The area consists of a recent alluvium layer 40 to 50 m thick overlying a Pleistocene 
formation that varies from 150 to 500 m in thickness. Underlying the Pleistocene material is a 
Miocene basement rock. Characteristic geological profiles also showing compressional wave 
velocities are shown in Fig. 2.3. Example soil profiles are shown on Fig. 2.4. The locations of the 


boreholes from which the profiles were constructed are shown on Fig. 2.5. 
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Fig. 2.1 Location of Lotung Large-Scale Seismic Test Site (LSST). 
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Five stages of laboratory testing programs were performed at the LSST during specific 
phases of the project to determine engineering properties of the soil. A summary of the tests 


performed is included here and the specific references for the test results are given below. 


1. 1984 Jong Shing Boring Services (JSBS) Laboratory Testing Program 


Index properties: 
Soil classification 
Grain size analyses 
Moisture contents 
Specific gravity 
Dry density and void ratio 
Atterberg limits 
Engineering properties: 
Direct shear tests 
Triaxial shear tests 
Unconfined compression tests 


2. 1987 National Taiwan University (NTU) Laboratory Testing Program 


Index properties: 
Grain size analyses 
Moisture content 
Specific gravity 
Dry density 
Atterberg limits 
Engineering properties: 
Uniaxial - load/unload and cyclic loading tests 
Triaxial - compression, extension and cyclic loading tests 
Resonant column 
Hydrostatic - load/unload and cyclic loading tests 
Compaction tests | 


3. 1987 University of California at Davis (UCD) Laboratory Testing Program 


Engineering properties: 
Triaxial shear loading 
Cyclic triaxial liquefaction testing 


The UCD phase of the laboratory testing was performed in conjunction with the installation 
of pore water pressure trausducers at the site. Pore water pressures were measured in the samples 


during these dynamic laboratory tests. 


4. 1989 NTU Laboratory Testing Program 


Engineering properties: 
One dimensional rebound 
Resonant column 
Cyclic triaxial liquefaction 
Permeability 
Cyclic triaxial modulus 


These additional tests were performed by NTU for the specific purpose of investigating the 
form of the shear modulus versus shearing strain and material damping versus shearing strain 


curves for undisturbed soil samples. 


5. 1990 UCD Laboratory Testing Program 


Engineering properties: 
Cyclic triaxial modulus 
Cyclic triaxial liquefaction 
Cyclic simple shear 
One dimensional rebound 
Permeability 


These tests by UCD were independent of the 1989 NTU tests and were performed to 
investigate discrepancies in modulus and damping data found from analysis of 1987 results for the 
LSST site. In addition, it provided additional data on the cyclic strength and liquefaction properties 
of the soil. Blowcount results from the standard penetration test (SPT) are shown on Fig. 2.6 for 
2 of the 3 drilling and sampling programs. The appearance of an occasional layer requiring an 
excess of 50 blows per foot indicates a gravelly soil zone that occurs in discrete lenses rather than 
as aconsistent layer. Results of cone penetration test (CPT) soundings are shown in Fig. 2.7 and 


2.8. Occasional spikes in tip resistance values confirm the presence of the gravelly lenses. 
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Fig. 2.3 Geological profiles for the Lotung site. 
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Fig. 2.4 Detailed soil profiles developed during the NTU field program. (EPRI, 1989) 
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Fig. 2.5 Location of boreholes and soundings at the Lotung site. (Anderson, 1993) 
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Fig. 2.6 SPT blowcounts recorded at Lotung test site. (Anderson, 1993) 
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Fig. 2.7 Cone tip resistance, friction ratio and estimated soil profile, CPT1. (Anderson, 1993) 
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Fig. 2.8 Cone tip resistance, friction ratio and estimated soil profile, CPT2. (Anderson, 1993) 
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2.4 Seismology 


A summary of the properties of measured temblors is given in Table 2.1. The epicenters 
for these events is shown on the map presented in Fig. 2.9. Events 4, 12, and 16 will be considered 
“large” events, with peak accelerations over 0.17 g. Events 12 and 16 were major events and have 
been discussed in detail (e.g. EPRI, 1989; Chang et al., 1991a, b; Anderson, 1993). Of these four 
large events, temblors 7 and 16 were deep focus events, and event 12 was an event occurring 
nearby and near the surface. Event 8, although a low acceleration event, is of special interest as it 


was an immediate aftershock of event 7. 


Examination of the Arias Intensity (Arias, 1970) for the events considered produces a 
different picture of the energy absorbed by the site for each earthquake. The Arias intensities for 
the Lotung temblors are given in Table 2.1, where it is seen that events 16, 12, 4, and 7 are the most 


energetic. The Arias Intensity, L,, is defined as 


BRE mw 
L, = 72) [4] dt, (2.1) 
0 


Where a is the acceleration time history. Whereas peak acceleration is simply the amplitude of the 
largest peak during the strong motion, Arias Intensity considers not only the amplitude of the 
acceleration but the duration and frequency content as well. Therefore Arias Intensity is a fuller 


measure of the disturbance at a given site. 


A reasonable assumption is that dynamic soil behavior is a function of excitation energy. 
This assumption has recently been supported by many careful experimental studies (e.g. Figueroa 
et al., 1994). The Arias Intensity is directly proportional to the energy input into the soil during 
earthquake excitation. In an attempt to normalize results to locally specific conditions, Arias 


Intensity and peak input acceleration will be used as a measure of event energy. 
2.5 Pore Pressure Generation 


Over the history of the Lotung experiment, more than 30 pore water pressure transducers 


were installed at the site (Shen et al., 1987). The location of these sensors in relation to the three- 
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arm surface array is shown in Fig. 2.10. The soil conditions at several pertinent locations are given 
in Table 2.2. As reported by Shen et al. (1987), most of the sensors ailed due to mechanical 
problems. However, several remained in operating condition and were triggered by events 12, 16, 
and 17. 

A typical pore pressure record is shown in Fig. 2.11 in relation to the acceleration time 
history. Information as to in situ pore pressure and increase for each sensor for events 12, 16, and 
17 are given in Tables 2.4 and 2.5, with the time histories from events 12 and 16 are given in 
Appendix C. 

Table 2.1 Properties of the recorded LSST series of earthquakes. 
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Note: The epicentral locations for Events LSST Nos. 13, 17 and 18 were not determined. 


Fig. 2.9 Locations of the epicenters for Lotung seismic events 1 through 18. (Tang et al., 1992) 
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Fig. 2.10 Location of the pore water pressure sensors at the Lotung test site. (EPRI, 1989) 
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Table 2.2 Soil conditions at the locations of pore water pressure sensors. (EPRI, 1989) 
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Fig. 2.11 Typical pore pressure record from event 16, LSST. Depth = 6.38 m; initial hydrostatic 
pressure = 72.6 kPa; increase in pore water pressure = 8.68 kPa. 
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Table 2.3 Peak pore water pressures, event 12. (Tang et al., 1992) 


Sensor 


No. 


PA-1 
PF-1 
Bra 2 
PE=5 
PN1-0 
PN1-1 
PN1-4 
PN2-1 
PA-Su 
PN3-1 
PN3-2 


Notes: 


(1) 
(2) 


(3) 
(4) 


Channel Depth, h 


No. (m) 

1 5.06 
5 335 
6 6.05 
9 12.00 
ney 316 
the 6.03 
aR) Shes) 
18 6.30 
22 Sisk 9, 
24 6.38 

29 E700 


u, (2) 
(kPa) 


54.60 
32.80 
95'170 
222.85 
Bb hgh 
61.30 


116.60 


Event LSST12 

(3) Aut 
(kPa) (kPa) 
59.90 5.30 
38.97 6.17 
5726] shah 
126.14 4.29 
42.12 2.40 
78.35 L7S 
65,53 Sinko 
736 36 11.24 
69.88 Lor S 
81.28 8.68 
123.41 6.81 


This table was modified from reference (2). 


u, = Hydrostatic pressure based on initia 


recordings. 


u, = Recorded peak pore water pressure. 
Au = u, - u, = Peak induced pore water pressure. 
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File 
Name 


E12COl. 
E12C05. 
E12C06. 
E12C0O9. 
EL2ZCi ae 
E12Ci23 
E12Ca 53 
E12C18. 
E12C21. 
E12C24. 
E12C297 


1 readings of the. 


Table 2.4 Peak pore water pressures, event 16. (Tang et al., 1992) 


Sensor Channel Depth, h Event LSST16 
No. No. ' (m) ae u? bul File 
(kPa) (kPa) (kPa) Name 
PF-8 HY 15.00 114.60 js Spe ey FAV TOs) E16C17.PWP 
PN2-1 18 6930 eume 66.90 83.60 16.70 E16C18 . PWP 
PA-3' 21 5.10 56.30 68.74 12.44  £16C21.PWP 
PN2-2' 23 8.00 97.20 106.61 = beet E16C23.PWP 
PN3-1 24 6.38 80.10 90.44 1O¥34 E16C24 . PWP 


Notes: (1) This table was modified from reference (2). 
(2) u, = Hydrostatic pressure base on initial readings of the 
recordings. 
(3) Se Recorded peak pore water pressure. 
(4) Au =u, - u, = Peak induced pore water pressure. 


at 


Table 2.5 Peak pore water pressures, event 17. (Tang et al., 1992) 


Sensor 
No. 


PF26 
PN2-1 
PA-3' 
PN2-2' 
PN3-1 


Notes: (1) 
(2) 


(3) 
(4) 


Channel 
No. 


ys 
18 
713 
23 
24 


Depth, h 
(m) 


15.00 
6.30 
3910 
8.00 
6.38 


u, ‘2? 
(kPa) 


114.60 


66.90 

56.30 
W720 

80.10 


Event LSST17 


“¢3) 
(kPa) 


128.73 


L2e47 
61.08 
LOT 92 
84.48 


Au’? 
(kPa) 


14.13 
5.87 
4.78 
4.72 
4.38 


This table was modified from reference (2). 
u, = Hydrostatic pressure based on initial readings of 
event SST16 recordings. 


u, = Recorded peak pore water pressure. 
Au = u, - u, = Peak induced pore water pressure. 
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File 
Name 


E17C17.PWP 
E17C18.PWP 
E17C21. PWP 
E17C23.PWP 
E17C24 .PWP 


2.6 Some Previous System Identification Studies of the Lotung Site. 
2.6.1 Zeghal et al. 


Zeghal et al., (1995) used cross-spectral analysis of the data set to calculate resonant 
frequencies up to the twelfth mode of vibration for the LSST site as a whole. The study showed 
consistent values for the six earthquakes used for the analysis. Elgamal et al. (1995) in a 
companion paper document a technique for evaluating the soil shear stress-strain histories directly 
from the downhole acceleration records using linear interpolation between downhole 
accelerometers. They use the results of the calculation to estimate the soil damping at different 


depths as a function of strain amplitude. 


Zeghal et al. (1995) estimated shear moduli and modal damping for a full suite of events, 
the results of which are included in Appendix D. Shear moduli were estimated for low-energy 
events (e.g. Events 5, 10, 13) to be 46 MPa for the 6 m depth, 63 MPa for the 11 m depth, and 69 
MPa for the 17 m depth. 


Zeghal et al. (1995) found much more scatter in their estimates of damping than for shear 
modulus, as expected from theory (this will be discussed in detail in Chapter 3), especially for low 
strain (< 3x10°%). Damping was also found to be at or above the commonly accepted values put 
forth by Seed and Idriss (1970). For Event 16, which generated an excess pore water pressure 
equal to about 25% of the effective overburden stress, the shear modulus reduction was about 14 


to 17%. 
2.6.2 Chang et al. 


Chang et al. have undertaken a series of studies of Lotung using parameter identification to 
evaluate the non-linear response of soils due to strong motion (Chang et al., 1991b, 1990, 1989). 
The first phase of the work concentrated on calculating a transfer function for the soil at several 
different points of the excitation history. The strong motion record was divided in three sections: 


(1) initial motion before strain levels high enough to degrade the soil was reached, (2) peak motion, 
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and (3) the coda or ring-down. The transfer functions were estimated from the ratios of smoothed 


output Fourier spectra and input spectra from various depths. 


The shear moduli showed a marked decrease from the initial low-level excitation to peak 
deformation. Depending on the depth used as input, the modulus reduction factor ranged from 0.60 
to 0.14. These values are for a magnitude 6.5 earthquake exhibiting a peak horizontal acceleration 
of 0.21 g. As acheck on reality, the S-wave amplitudes for the initial segment were in very good 
agreement with the actual measured values. A troublesome point is the large discrepancy in 
modulus reduction factor for the two horizontal components. The expectation is that they would 


be virtually identical, since soil degradation is usually thought of as a scalar quantity. 


The results from the forward propagation non-linear DESRA-2 (Lee and Finn, 1985) 
analysis show good agreement between actual and calculated displacement for frequencies up to 
about 6 Hz. There was also good correlation for phase angle for frequencies up to about 3 Hz. The 
shape of the equivalent damping ratio curve back-calculated from the non-linear analysis is 
different from that normally expected. Rather than the expected hyperbolic curve, the field curve 

S-shaped. This uncommon S-shaped curve is the same shape as that calculated by Abdel- 
Ghaffar and Scott (1979) from actual field data. 


Further work on the shear modulus reduction curve based on actual large strain 
measurement is reported by Chang et al. (1991b), using shearing strains caiculated by SHAKE. 
The results show that for small strains (surface acceleration < 0.03 g), the results from geophysical 


methods, resonant column, and system identification are in good agreement. For intermediate 
: -3 -2 . 
strains (5x10 to 2x10 °) the back calculated values for modulus reduction are up to twenty percent 


lower than the resonant column measurements. For higher strains of 3x10 to 2x10", the results 
from cyclic triaxial tests are in fair agreement with the field values for shear modulus, with 


moderated scatter for the laboratory data. 
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CHAPTER 3 - PROCEDURES FOR SYSTEM IDENTIFICATION 


eva System Identification 
3.1.1 Parametric Modeling 


The goal of system identification is to model a system in a manner that provides needed 
mechanical information about that system. The most common techniques have evolved from 
electrical and mechanical engineering, and involve solving the inverse problem for the system 
transfer function. Each method has limitations; in the words of G. E. P. Box, "All models are 


incorrect, but some are more useful than others." 


The process of inversion allows the estimation of the system response function (filter) if the 
input and output signals are known. A simple model for characterizing a system is as a parametric 
relationship between system input and output. Such a model, referred to as an autoregressive- 


moving average (ARMA) model, is based on discrete time series analysis: 


nb na 
Oe pe Go eo te ct Og tt DX t ap D Xp Da apy ek (3-4) 
j=0 k=1 


where y; is the actual output data sequence, x, 1s the input sequence (assume white noise for simple 


spectral estimation), na and nb are the AR and MA orders, respectively, and the subscript is the 
time step counter. The output is seen as a combination of the input history acted upon by the "b" 
coefficients plus the past outputs acted upon by the "a" coefficients. The input series, involving 
the "b" coefficients, is a causal moving average (MA) process (convolutional). The series 
involving weighted past output values ("a" coefficients) is a noncausal autoregressive (AR) 
process. The lengths of the AR and MA processes (model order) must be explicitly chosen so that 


the model best represents the process. 
Applying the shifting theorem to Eq. 3.1. yields the Fourier transform (Bracewell, 1978) 
2i@ 


Vell tale. ia,es wee eeeNM (EN been hse) t..) (3.2) 
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where i is ./—1 and @ is circular frequency. Applying the Z-transform (Bracewell, 1978), where 


z = 2? , and rearranging, yields the frequency domain transfer function H, 
1 2 
Y 0, + Oe ee 
De = (3.3) 
X 1 2 
a) mens wl LYS ns cow 


The ARMA model is very powerful in that it can easily model sharp drops, sharp peaks, 
and smooth spectral behavior. It is also the most parsimonious estimator (Robinson, 1982), 
describing a complex process with very few parameters calculated from a small length of data. 
Parametric modeling avoids many of the difficulties inherent in the traditional Pouce methods, 
discussions of which can be found in many books and journals (e.g., Glaser, 1993; Johansson, 
1993; Pandit, 1991). Extensions of this model, e.g., ARMAX, ARX, Box-Jenkins, allow input, 
system, and output noise to be expressly modeled (Ljung, 1987). In particular, the ARX model 


includes the effect of uncertainties and noise as a white noise term. 


The ARMA model has special significance since it can be derived directly from the 
differential equation of motion for an N-degree-of-freedom (DOF) system, with the damping ratio 
and resonant frequency as the model parameters (e.g., Gersch and Luo, 1970). A 2n-2n ARMA 
model is therefore a valid model for a layered soil system, or soil-structure interaction problem. 
The damping ratio and resonant frequency of the N-degree-of-freedom oscillators are contained in 
the 2n AR parameters. Phase relations are preserved in the MA parameters. The modal frequencies 
€;, percent of critical damping ;, (Ghanem et al., 1991) and power participation factor p; (Pandit, 
1991; Safak, 1988) are calculated from the system poles and residues found from partial-fraction 


expansion of Eq. 3.3. The modal parameters are defined as 


_ Aree 


= 3.4 
QO; AG (3.4) 
}. 
E> a (3.5) 
(hi + 8 
pj = —r,conj(z,) —z ;Conj(r ;) (3.6) 
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where A = Arg(z pod = —(0.5)In|z is , Z is the pole for mode j, r; is the residue for mode j, and 


At is the digitization rate. 
3.1.2 Adaptive (recursive) model estimation 


Traditional methods of system estimation, both parametric and non-parametric, are strictly 
valid only for stationary data. A stationary signal is one whose statistics do not change with time. 
The commonly invoked, loose definition of stationarity requires that the variance of the signal be 
constant over any and all time windows. Inherent in the Fourier transformation of a time series to 
the frequency domain is the averaging of the signal components over the sampling period T. A 
piece of time is frozen over this period and the assumption made that all time before and after is 
the same, i.e. repeated forever. The energies present at each component frequency are integrated 


over the entire time period T. 


The difficulty with non-stationary signals is that these energies are changing during this 
period. If the frequencies present are changing over this time window, the resulting estimation, 
regardless of method used, will be a smeared average as if all the frequencies with a given energy 
were active throughout the entire period. For weakly non-stationary processes, the effect over a 
small time period is unimportant. If needed, the signal can be cut into relatively stationary sections 


and spectra found using methods specially designed for short data segments, i.e. Burg's method. 


The field of adaptive filtering was formed to model non-stationary processes. As the 
Statistics of the signal change through time, the filter "adapts" to the changing variance with new 
parameters that reflect the structure of the system at that point. The predicted value for the next 


time step can be compared with the actual value, and the difference (referred to as innovations) 


sh a 
(y,- 91) (3.7) 


where y, is actual output at time t, and j, is the prediction of output at time ¢ made at time f-J, will 


give a measure of how well the filter is doing its job. The term "innovations" is used because this 


information is new information that can not be predicted by the model at this particular step. 
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Autoregressive parameters can be sequentially estimated so that the parameters are 
adaptive to the changing nature of the process (Marple and Lawrence, 1987). The AR parameters 
are updated after each data point, tracking slowly non-stationary signals. A forgetting factor, 
commonly a damped negative exponential, is used so that older data carries less and less weight. 
A frequency domain estimation can be made at any time step by evaluating the AR parameters 


around the unit circle, giving the spectral representation of the behavior of the process at that time. 


The most popular direct adaptive filter, or process model, is the so-called Kalman filter 
(Kalman, 1960; Kalman and Bucy, 1961). Sorenson (1970) points out that the Kalman approach 
is a direct descendant of Gauss's least squares, except now neither the signal nor the noise model 
must be stationary — the state may change from sample point to sample point. Nau and Oliver 
(1979) state that the Kalman filter is based on a dynamic AR model defined by "two concurrent 


random equations of motion": 


x, = H’,_1®,+a, (3.8) 


the AR(p) equations of motion, and the "motions" of the parameters, 


®, = 9,_, +8, (3.9) 
where p = number of prior observations utilized, 

Hs = vector of p prior data observations KS A DS 

, = vector of p AR parameters, 

at = Gaussian white noise with 0 mean and variance o” 

by = Gaussian white noise with 0 mean and covariance matrix Q. 


Equation 9 estimates a value of ©, comprised of p previous parameters, through a random 


walk equation. The estimate uses the weighted p previous data points, and yields a new 


observation x; when added to a new noise value. The least squares solution solves the equations 


so that the innovations (Eq. 3.7) — new, dynamic information that cannot be predicted — are 


minimized in a least squares sense each time step. 
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The theory behind the Kalman filter can be manipulated to yield the system parameters for 
the case where there is no a priori information about the noise, and even when there is no 
information about the input signal. The so-called extended Kalman filter has been very 
successfully applied to non-stationary (and non-linear) estimation problems (Ljung, 1979; Astrom 
and Eykhoff, 1971). The manner of application is actually straight-forward. The Kalman model 
is constantly updating its estimation of the dynamic process by examining the innovations. The 
dynamics can be due to a changing input or noise process, or it can be due to the system itself 
changing. The effect is a linearization between single time steps, but if the system is changing 
slowly compared to the time step used, the linearization is "invisible" and the non-linear behavior 


is well modeled. 
3.2. Application of Parametric Models 


Previous work sponsored by NIST (Glaser, 1995) demonstrated that the MATLAB 
software package (MathWorks, 1993) is ideally suited for manipulation, processing, and 
presentation of earthquake data. MATLAB is a matrix-based system which evolved from the 
LINPAC and EISPAC libraries commonly used for mainframe FORTRAN numerical analysis. 


Complex numerical problems can be speedily solved without programming in the traditional sense. 


For this study, the standard routines contained in the MATLAB System Identification 
Toolbox (Ljung, 1993, 1987) allowed SI to be used as a tool accessible to the geotechnical 
engineer. Virtually every approach and algorithm encountered in the literature by the author could 
be duplicated rapidly and accurately. When run on an Silicon Graphics Indy workstation, all 


aspects of the analysis were quick enough to be interactive. 


Analysis begins by determining if the event can be modeled as stationary segments. A 
recursive segmentation scheme, which attempted to break the data into segments with a chosen 
maximum variance (Ljung, 1987) is used. However, it is not possible to determine the “correct” 
variance a priori. In these cases a more direct method is used — the output simulated by the 


calculated system has to accurately model the actual measured output. 


ID, 


3.2.1 Choice of Model Order, and Validation 


When possible the input-output data record is broken into segments based on a mechanistic 
understanding of the seismic event and soil behavior gained from study of the pore pressure 
behavior. Initially it is assumed that the various segments are basically stationary. If the stationary 
model can not accurately and parsimoniously simulate the segment output, a non-stationary 
recursive model is used. In addition, the appropriateness of the model is checked by insuring 99 
percent confidence in both the whiteness of the residual autocorrelation function and the cross- 


correlation function between the input and output residuals (Bohlin, 1987). 


The stationary SI algorithm uses a least squares estimation for the ARX model. It is 
necessary to estimate the number of parameters to be calculated, which is essentially estimating 
the degree of freedom of the soil system. There is no obvious answer to the degrees of freedom of 
the system, so several verification techniques are employed to insure that a proper model order is 
estimated. Model order is increased in 2n-2n steps corresponding to an additional degree of 
freedom each. The simulated output of the model is then compared with the actual output for 
congruence, and the fewest parameters needed to accurately characterize the system was chosen as 
the system model order. Examination of the pole and zero plot insures that excessive, overlapping 
parameters are not included (Ljung, 1987). If the segments proved non-stationary, they are 
analysed using a tecursive Kalman filter technique that expressly accounts for non-stationarity 


(Ljung, 1987). The trade-off relative to the stationary ARMA algorithm is reduced certainty. 


Each input/output pair was initially submitted to an overall algonthm that calculates the 
loss function (normalized sum of squared prediction errors, Eq. 3.7) versus model order for a suite 
of pre-selected model orders. A typical loss function plot is shown in Fig. 3.1, for Event 16, DHA 
0-11 m. As is obvious from this plot, the waveforms only carry information from the first 3-4 
modes (12-16 parameters) - there is very little improvement in estimation quality for more 
parameters. The results of the loss function calculation were compared with the model orders 
chosen by other criteria, in particular the Akaike information theoretic criterion (AIC) and 
Rissanen's minimum description length (MDL) tests, which gave the same results. The AIC looks 


for a mode] that minimizes the 'information distance’ to the actual system. The MDL searches for 
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a mode] that minimizes the length of the description of the data, i.e. the most efficient model 


(Ljung, 1987). 


In most cases, the fits were excellent for a smal] number of parameters. An example is 
shown in Fig. 3.2 which shows the model and actual values of output for DHA Event 16, 0-11 m, 
NS component for a 2-DOF model (8 parameters). For this calculation, the rms error 
(unnormalized loss function) varied from 1.037 for the 2-DOF model to 1.018 for a 10-DOF 
model, indicating that no new information is presented by the more complicated models. In 
general, the match between actual and modeled outputs for all waveforms investigated were not 
fundamentally improved by using up to 40 parameters (10-DOF). Given that the data does not 
have an infinite signal-to-noise ratio, i.e. there is noise present from many sources including 
quantization, there is only a limited amount of information that can be taken from the data (e.g. 
Pierce, 1980; Shannon, 1949; Shannon and Weaver, 1949). Experience has shown that any 
estimate above the second mode is tenuous at best, although numerically we are only limited by 


computational power as to how many modes we want to calculate. 


The decision-making process as to model choice will now be shown in detail. Figures 3.3 
to 3.13 document the process of model selection for a rapidly converging choice. Figure 3.3 isa 
plot showing the input and output data set for DHA Event 12, 6 m to surface interval. From the 
similarity of the input and output time series we would expect that a relatively uncomplicated filter 
would describe the system. Figure 3.4 is the loss function plot for that interval. The steep drop in 
loss function between 4 and 8 parameters is another indication that a simple system description is 
expected. Figures 3.5, 3.6, and 3.7 show the comparison of model vs. actual output for 1-DOF, 3- 
DOF and 5-DOF systems, respectively. Examination of the curve fits for the various models 
indicates a slight improvement in fit between actual and estimated interval output with increasing 
model order, from 0.3459 to 0.3056. In fact, the change from the 3-DOF to the 5-DOF system only 
improves the match by 0.0055. 


The residual analysis plots for the same models, shown in Figs. 3.8 - 3.10, also demonstrate 
that no additional statistically significant information is being extracted beyond the 3-DOF system. 


Comparison of Figures 3.11-12, the pole-zero plots for the 3-DOF and 5-DOF models, show that 
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Figure 3.1 Loss function versus number of modal parameters used to model Event16 N-S, DHA, 
O - 11 meters depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 1.037 
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Figure 3.2 Comparison of actual and modeled output for Event16 N-S, DHA, 0 - 11 meters depth 
interval. 
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Solid: Input signal, Dashed: Output signal 
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Figure 3.3 Comparison of actual signal input and output for Event16 N-S, DHA, 0 - 6 meters depth 
interval. 
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Figure 3.4 Loss function versus number of modal parameters used to model Event12 N-S, DHA, 
0 - 6 meters depth interval. 
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Figure 3.5 Actual versus 1 - DOF model output for Event12 N-S, DHA, 0 - 6 meters depth 
interval. 
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Figure 3.6 Actual versus 3 - DOF model output for Event12 N-S, DHA, 0 - 6 meters depth 


interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.3138 
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Figure 3.7 Actual versus 5 - DOF model output for Event12 N-S, DHA, 0 - 6 meters depth 
interval. 
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Correlation function of residuals. Output # 1 
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Figure 3.8 Residual analysis plot for 1-DOF model, DHA, event 12, 0-6m, n-s 
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Correlation function of residuals. Output # 1 
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Figure 3.9 Residual analysis plot for 3-DOF model, DHA, event 12, 0-6m, n-s 
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Correlation function of residuals. Output # 1 
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Figure 3.10 Residual analysis plot for S-DOF model, DHA, event 12, 0-6m, n-s 
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Figure 3.11 Pole - zero plot for 3 - DOF model of Event12 N-S, DHA, 0 - 6 meters depth interval. 
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Imaginary Component 


Real Component 


Figure 3.12 Pole - zero plot for 5 - DOF model of Event12 N-S, DHA, 0 - 6 meters depth interval. 


43 


the poles and zeros begin overlapping for the 5-DOF case, due to no “new information” available 
to be extracted, indicating the model order is too high or unnecessarily complicated for the interval 
being studied. Thus a 3-DOF model was the choice to represent this interval and event. Table 3.1 
summarizes the calculated modal frequencies and damping ratios for the models examined in this 


analysis. 


In some cases the models generated did not adequately represent the system for a very small 
model order. Figure 3.14 is the loss function plot for DHA, event 4, 6-11 m depth, east-west 
component showing a slower rate of convergence with increasing model order compared to the 
previous example. The fit for the 1-DOF model is 0.6347. The curve fits initially become poorer 
as the model order increases to a maximum of 0.702 for the 3-DOF system. From a 4-DOF to a 5- 
DOF model the fit improves again from 0.6852 to 0.6185. Figures 3.15-16 are the residual analysis 
plots for the 3-DOF and 5-DOF systems respectively, and indicate an improvement in information 
extraction is occurring for the 5-DOF model over the 3-DOF model. A 6-DOF model, however, 
shows a decrease in the fit of the actual to predicted values to 0.6317 without additional 
information utilization. This can be demonstrated by comparing Figure 3.17, the residual analysis 
plot for the 6-DOF model, to Figure 3.16. For this interval a 5-DOF model was chosen to represent 


the system for this event. 


Some intervals could not be well matched using a standard ARX model. These intervals 
were successfully fitted using either a recursive RARX model. Recursive techniques allow the 
system parameters to change every time step if necessary to produce a suitable match of model to 
actual system output. These algorithms are thus ideal for non-stationary or difficult to match 
signals. Figure 3.18 is a plot showing the best fit of an ARX model to actual output for DHA, Event 
4, 47-17 m depth. In this case the “best-fit” ARX model is a 12-DOF or 48 parameter model! 
Figure 3.19 shows the increased performance in fit for a recursive model, in this case an 3-DOF 


RARX model. Note that this implies that the system is changing - softening - through time. 


Table 3.1 Calculated modal frequencies and damping ratios for models calculated for Event 12 N- 
S, DHA, 0 - 6 meters depth interval. 
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Figure 3.14 Loss function versus number of parameters, Event 4, DHA, E-W, 6 - 11 m depth. 
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Correlation function of residuals. Output # 1 
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Figure 3.15 Residual analysis plot for 3-DOF model, Event 4, DHA, E-W, 6 - 11 m depth. 
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Correlation function of residuals. Output # 1 
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Figure 3.16 Residual analysis plot for 5-DOF model, Event 4, DHA, E-W, 6 - 11 m depth. 
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Correlation function of residuals. Output # 1 
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Figure 3.17 Residual analysis plot for 6-DOF model, Event 4, DHA, E-W, 6 - 11 m depth. 
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Velocity (cm/sec) 


Figure 3.18 


Event 4, DHA, 47 m to 17 m, arx best-fit [12 12 2] 


——-" actual 
‘—-- estimated 


10 ihe) 20 25 30 33 40 
Time (sec) 


“Best Fit” ARX model for Event 4, DHA, 47 m to 17 m depth interval, 6-DOF. 
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Figure 3.19 Recursive RARX model for Event 4, DHA, 47 m to 17 m depth, 3-DOF. 
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3.3 Meaning of the Estimated Damping Values 


A question must be raised about the meaning of the damping values calculated in this, and 
other, studies — at this point the geotechnical engineer does not really know what the “correct” 
answer should be! A common practice is to hold up “damping” values that work well in SHAK 
(Schnabel et al., 1972) as the “real” answer, but these values are also a numerical “nicety” to make 
a given model work. The parameters calculated in this study represent the weights that can transfer 
input into output through convolution. These parameters can be represented in the frequency 
domain as poles and zeros, and in turn as resonant frequency and damping. The detailed 
congruence between actual system outputs and modeled outputs will be discussed in this chapter, 
and are presented for all studied events in Appendix E. The values calculated, reasonable 
compared to other modal studies (e.g., Iemura et al., 1990) and many laboratory studies (Vucetic 
and Dobry, 1991), are better seen as “effective” damping coefficients which capture all forms of 


frequency-dependent mechanical losses, rather than an ‘intrinsic’ material property. 


A problem common to all estimation techniques is that the damping is being modeled as 
viscous damping of a lumped-mass system because it is the most amenable to calculation (Pandit, 
1991). In actuality “damping” is measuring cumulative energy dissipation which will include 
viscous damping, plastic. deformation, friction and attendant heat generation (Kramer, 1995). Even 
with a perfect technique, the damping values estimated by the model will not truly represent what 


is physically occurring in the soil 


Damping estimates are often very sensitive to subtle changes in the modeling of the 
system, especially for recursive estimates which have very different system parameters changing 
rapidly. In this case a limited number of data points enter directly into the calculation and variance 
is inversely proportional to the square root of the number of data points. The physical 
interpretation of the instantaneous mechanical values is also not immediately clear. The ARMA 
parameters recursively c:..culated at any given time define the filter needed to transform that time 


step of input data into the » >xt output value. The mechanical characteristics are mathematically 
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extracted from the filter based on a limiting set of assumptions, and an “instantaneous frequency” 
or damping has no physical meaning. The ARX estimates are made with more input data points 


than the Kalman estimates, but the Kalman assumptions better model the changing system. 


Damping is inherently harder than natural frequency to identify by any method or 
calculation. A study was done by Gersch (1974) in order to determine the greatest degree of 
accuracy with which a proper order ARMA model can estimate the damping ratio and natural 
frequency of a structure, using the Maximum Likelihood method. He notes that as the number of 
data points (N) becomes large, the estimates approach the actual, and the model errors approach 
the Cramer-Rao lower bound of variance. For both parameters, the coefficient of variation is 
inversely proportional to N and length of the period sampled, and relatively insensitive to noise. 
For one thousand data points, the coefficient of variation was less than 0.01 for natural frequency, 


but greater than 0.2 for damping ratio; this is an order of magnitude difference. In summary: 


* error in damping estimates will be one to two orders of magnitude larger than for natural 
frequency. 


* error for both factors are inversely proportional to the number of data points included in 
the analysis and the sampling period 


* error for both factors are insensitive to additive noise when a MLE procedure is used and 
the number of modes of the system being modeled 


3.4 The Pseudo-Stiffness 
As the pore pressure ratio increases, the effective stress acting among the sand grains 
decreases, and the shear modulus (G) decreases. Since it has been impossible to accurately 


measure stress and strain inside the soil body, the ratio of applied force to displacement (stiffness) 


is used to characterize material behavior. 


The here-proposed pseudo-stiffness is calculated by the following algorithm. 


1) All horizontal motions from the surface (s) accelerometers are accounted for by 
taking the Euclidean norm of the north-south and east-west surface acceleration 


rd. =a : a 
récoras — ans = a, (N-S) + 4s (E-W) : 
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2) All horizontal motions from the buried (b) accelerometers are accounted for by 
taking the Euclidean norm of the north-south and east-west buried acceleration 


, roy ey 2 


at The relative acceleration of the surface to buried transducer is determined from the 
step-by-step difference between the surface and buried horizontal accelerations — 


a = @ ==! (Cl 
reli) Sa hi) b, hy) 


4) The relative accelerations are converted into forces by multiplying the differential 
accelerations acting at the top and bottom of the layer of interest by the soil column 
mass, which should change relatively little during the softening process (e.g. 


Udwadia, 1985) = re = Graliy : cy. * A, 


5 The running total of the amount of force is computed for the system during each 
time step — sum(t)=sum(t-1)+{force(t)-force(t-1)} — yielding the cumulative 
horizontal unit force. 


6) Steps I, 2, 3, and 5 are applied to the displacement records to yield the cumulative 
horizontal displacement. | 


7) The ratio of the resulting cumulative horizontal unit force to displacement are two 
independent variables that combine to give a pseudo-stiffness. 


The use of a single transducer to provide acceleration and displacement measurements is 
not problematic since the derived displacement (Brady et al., 1989) is generally assumed to be 
identical to what would be recorded by an independent displacement transducer placed next to the 
accelerometer. The calculations do not distinguish between active (external, inertial) and 
restorative (internal, stiffness) forces acting on the soil. The pseudo-stiffness is a convenient 
representation of the system behavior, combining load, time, and displacement information in a 


simple curve. 
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CHAPTER 4 PARAMETRIC MODELING OF LOTUNG DATA 


4.1 Introduction 


As discussed in Chapter 2, eight earthquakes of varying characteristics were chosen for this 
study. The selected temblors represent two classes - large and small. The estimated mechanical 
parameters for the small (low energy/acceleration) events can be seen as a linear baseline against 
which to judge whether the larger events exhibited nonlinear behavior. This chapter will begin 
with a general discussion of the results of the study, followed by general discussions on modal 
frequency, damping, and soil-structure interaction for the lotung site as a whole. The two similar 


small events (events 7 and 8) and the large events (events 4 and 16) will then be discussed in detail. 


All pertinent intervals for the Lotung events were modeled during this study. The 
comparisons of all actual interval outputs to estimated system outputs are presented in Appendix 
E. An example of the congruence between estimated and actual output is given in Fig. 4.1. This 


result is for the top six meters for the free-field DHB array due to Event 4, using a 2-DOF model. 
4.2. Discussion of Results 


A summary of the estimated fundamental frequencies for the earthquakes analyzed for this 
report is given in Table 4.1. The corresponding estimates of system damping are given in Table 
4.2. As demonstrated in Tables 4.1 and 4.2, when subjected to greater seismic energy the modal 
frequencies do systematically decrease and damping ratios increase. In fact, the results strongly 
support the concept of threshold strain (Dobry, 1973). The modal frequencies and damping ratios 
are very close for small Events 3, 8, 9, and 10. These events all have values of Arias Intensity 
under 200 cm/sec and maximum acceleration values of 0.07 g and less. The higher energy events 
including Events 7, 4, 12 and 16, with Arias Intensities between 650 and 4000 cm/sec and 
maximum accelerations ranging from 0.17 to almost 0.5 g exhibit substantial and consistent 


decreases in fundamental frequency. 


a. 


velocity (cm/s) 


Event4, DHB, E — W, surface to6m 
40 


Solid: Model output 


30 Dashed: Measured output 


RMS Error: 0.9253 


20 
10 
0 
-10 
—20 
0 5 10 15 20 25 30 oo 40 
Time (sec) 


Figure 4.1 An example of “goodness-of-fit” between actual and modeled system output. 
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Consider the results from the "free field" hole, DHB. The fundamental frequency for the 
0-6 m interval exhibits a decrease from 4.6-5.2 Hz for Events 8, 9, 10 and 3 (Arias Intensity ranging 
from 36-192) to 2.6-3.9 Hz for Events 7, 4, 12 and 16 (Arias Intensities 652-3995.) This result is 
consistent with the calculations from deeper intervals. The 6-11 m segment shows a decrease from 
an average of 3.3 Hz for the moderate events to 1.4-2.4 Hz for the larger events. The trend 
continues to the 11-17 m interval where the fundamental frequency decreases from 2.6 Hz to 1.6 
Hz. For DHB the only deep interval information is the result from Event 4, which also exhibits a 


decrease in frequency from 1.4 Hz to 0.9 Hz for the 17-47 m segment. 


Figure 4.2 graphically compares estimates of fundamental modal frequency to peak 
acceleration for all events, while Fig. 4.3 graphically compares estimates of system damping to 
peak acceleration for all events. The decrease in fundamental frequency and increase in damping 
with increase in peak acceleration is readily evident. Close examination shows that the change in 
material properties is greater for the near-surface layers than for the deep one. It is expected that 


greater confinement of deeper soils would serve to smooth out changes in system behavior. 


Figure 4.4 graphically compares estimates of fundamental modal frequency to Arias 
intensity input into the soil layers, while Fig. 4.5 graphically compares estimates of system 
damping to Arias intensity input into the soil layers. Surprisingly there is less correlation between 
system behavior and energy input into soil layers as measured by Arias Intensity. This leads to a 
conclusion that soil softening is more a function of peak acceleration (which causes peak strain) 
than total input energy (that is proportional to cumulative strains). This conclusion implies a 


minimum peak acceleration is needed to excite the soil past its threshold strain. 


The estimates for system damping presented in this study are greater than values 
traditionally accepted by geotechnical engineers (see Appendix D). Estimated values range from 
less than 2% to almost 50%. As discussed in Chapter 3, these are the values that map system input 
into system output. A quick review of the estimation results in Appendix E show how well the 
estimation process was able to capture the essence of the soil systems in question. Damping values 
reported in Table 4.2 are generally higher than estimates of damping by Zeghal et al. (1995) for 
Events 7, 12, and 16 (see Figs. D.9 through D.18) 
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A summary of second mode frequency estimates is given in Table 4.3, and damping 
estimates in Table 4.4. The damping is in general higher than for the fundamental modes, although 
the there is a fair amount of scatter due to the estimation process. If a system needed more than 2- 
DOF to model the soil response well, the second mode is often a “phantom” mode that does not 
represent strong system behavior. Therefore, the participation factor must be examined to 


determine whether the second and higher modes are significant. 
4.2.1 Soil-Structure Interaction 


The basic understanding of soil-structure interaction (SSI) predicts two effects of a 
structure sitting on a non-rigid foundation: (1) the fundamental frequency of the soil/structure 
system would be lower than that of the structure on a rigid base and (2) the damping ratio of the 
soil structure system will be larger than that of the structure on a rigid base (Kramer, 1995). These 
effects are in terms of the building rather than for the soil itself. Rollins and Seed (1990) discuss 
the varicus changes to soil response possible but come to no firm conclusions. The same is true 


for the extensive study on SSI for the Lotung site (Tang, 1987). 


Damping and modal frequency results “or DHA are given in Tables 4.1, 4.2, 4.3, and 4.4. 
The res’ *s are summarized graphically in Figs. 4.2 through 4.5 in cor <arison to the free field 
results. 1. 1s interesting to note that the decrease in fundamental frequency is much less pronounced 
for the results from DHA, the array under the model containment structure. For the 0-6 m interval 
the fundamental resonant frequency shifts only from 3.1 Hz for Event 3 (maximum acceleration 
0.01g) to 2.2-2.4 Hz for Events 4, 12 and 16 (maximum accelerations of 0.17-0.49 g). Note also 
that the fundamental frequencies for the same events in the 0-6 m interval are lower under the 
structure compared to the free field values. T’ ‘s difference is no longer evident by the 6-11 m 
level. In fact, examination of the data for the 6-11 m, 11-17 m, 17-47 : and 11-47 m intervals 
indicate hat modal frequencies are equivalent at DHA and DHB below 6 m. Other authors have 


predicted that the soil structure interaction wo id be evident to a depth of 17 m (Tang, 1987). 


Examination of pseudo-stiffness provides another method of examining soil-structure 


interaction. While detailed analysis of this eff«ct will follow in the sections on the individual 
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earthquakes, the results showed a clear effect for the top interval (6 m to the surface) and in rare 
occasions for the second interval (11 m to 6 m). For Events 4, 7, 8, and 16 the pseudo-stiffness is 
smaller at the structure than in the free-field. For event 12, the surface layer of soil appears to have 


a greater pseudo-stiffness than in the free-field. 
4.3 Event 4 


Event 4 subjected the Lotung site to the highest peak acceleration of any recorded event in 
the LSST suite (0.49 g), although the strong motion for Event 4 was very brief (about 1.6 seconds) 
as opposed to event 16 (about 12 seconds). This lead to a much lower Arias Intensity (1845 m/s 
vs. 3995 m/s) and a slightly smaller magnitude (M, = 6.5 vs. 7.0). Because of the high peak 
acceleration and other properties (Elgamel, 1995) researchers thought that the pore water pressure 
generated during this event must be at least as great as that recorded during event 16. Since there 
is no pore water pressure record for Event 4, other methods must be found to estimate system 


behavior. 


The Event 4 surface layer (6 m to the surface) signal from DHA proved difficult to model 
using stationary techniques, indicating some degree of non-linear behavior that can be captured 
using recursive techniques. It was found that the system could be well modeled as a 3-DOF ((6 6 
1]) system using a forgetting-factor of 0.83 (Ljung, 1987). Figure 4.6 shows the congruency 
between the actual and modeled system, while the evolution of the three modal frequencies (1 Hz, 
3.5 Hz, and 7.3 Hz) is given in Fig. 4.7. It is seen that the system behaves very erratically for the 
first 10 seconds, or the period of strong shaking. It is evident from the velocity history shown in 
Fig. 4.6 that the signal becomes more regular and lower in frequency content at this point. The 
change is reminiscent to that of a soil undergoing liquefaction, such as the time histories from the 
Wildlife site, Superstition Hills, earthquake (e.g. Glaser, 1996, 1995). The evolution of the Event 


4 surface layer system damping could be estimated, and is given in Fig. 4.8. 
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Figure 4.9 Event4, DHB vs. Event8, DHB; 6 m to Surface 
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Figure 4.10 Event4, DHB vs. Event8, DHB; 11 mto6 m 
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Figure 4.11 Event4, DHB vs. Event8, DHB; 47 to 11m 
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4.3.1 Event 4 Pseudo-Stiffness and Soil-Structure Interaction 


One method is to compare the pseudo-stiffness behavior of the two events to determine if 
there was comparable softening taking place, in particular a comparison with a very small event 
such as Event 8. Figure 4.9 shows a comparison of the pseudo-stiffness for the deep interval, 47 
m to 11 m (due to the confinement it is expected that this interval would show no loosening, as 
shown in Sect. 4.2). It is seen in detail that initially the pseudo-stiffness is identical for both Event 
8 and 4, with slight reduction in pseudo-stiffness after the soil absorbed most of the Event 4 
energy. The same observation holds true for the 11 m to 6 m interval shown in Fig. 4.10 and the 


6 m to surface interval shown in Fig. 4.11. 


However, when Event 4 pseudo-stiffness is compared to Event 16 pseudo-stiffness, the 
values are the same. Figures 4.12, 4.13, and 4.14 show a comparison in pseudo-stiffness between 
events 4 and 16 for the 47 m - 17 m, 11 m- 6 m, and 6 m to surface intervals, respectively. It is 
readily seen that the slopes shown (the pseudo-stiffness) are the same, although there was more 


total motion during Event 16 (AI= 4000 m/s) than Event 4 (Al= 1850 m/s). 


Comparison of the pseudo-stiffness results from the free-field array (DHB) and the array 
adjacent to the containment vessel (DHA) show effects of soil-structure interaction. Figure 4.15 
shows the pseudo-stiffness curve for both arrays, where it is seen that the slope is identical for the 
two. When the behavior for the 11 m to 6 m interval is examined in Fig. 4.16, there is evidence 
that the structure acted to smooth and constrain the soil behavior, although the soil ends up 
behaving the same in both locations. For the top interval, 6 m to the surface shown in Fig, 4.17, 


the soil adjacent to the structure shows significantly greater softening than that in the free-field. 


This can be compared with values given in Table 4.1 for fundamental frequency. There is 
only a significant difference in fundamental mode for the top, 6 m to surface, interval - about 2.2 
Hz at DHA and 2.8 Hz in the free-field. The damping estimates, from Table 4.2, are too erratic to 


yield any firm conclusions. 
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Figure 4.12 Event4, DHB vs. Eventi6, DHB; 11 mto6m 
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Figure 4.13 Event4, DHB vs. Event16, DHB; 11 mto6m 
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Figure 4.14 Event4, DHB vs. Event16, DHB; 47 m to 17 m 
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Figure 4.15 Event4 free—field vs. Event4 SSI; 6 m to surface 
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Figure 4.16 Event4 free—field vs. Event4 SSI; 11 mto 6m 
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Figure 4.17 Event4 free-field vs. Event4 SSI; 47 m to 17m 
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4.4 Events 7 and 8 


Event 8 (M,=6.2) is considered to be an aftershock of Event 7 (M,=6.2) (see Section 2.4), 
So comparison of the results from these two events should cancel out any effects associated with 
hypocentral location and source parameters. Although both events had very similar magnitudes, 
peak acceleration for Event 7 was much greater than for Event 8 (0.21 g vs. 0.03 g) as was the Arias 


Intensity (652 m/s vs. 139 m/s). 


However, the models for the low-energy Event 8 was able to capture the behavior of the various 
intervals well enough to accurately estimate the behavior of the Event 7 shaking. Figure 4.18 
compares the actual Event 7 surface layer output to the output estimated from the appropriate Event 
8 filter. The input to the model was the N-S 6 m velocity time history. Further evidence as to the 
model capturing the essence of the layer is given in Fig. 4.19 which is a comparison of the Event 


7 E-W output to the Event 7 output modeled through the Event 8 N-S filter. 


Figure 4.20 shows the results of modeling Event 8 N-S output through the Event 7 N-S 
filter. The filter captures the frequency information, with only a very slight (2 time step) phase lag, 
but tends to overestimate the early, high, amplitudes by about 10%. These results are consistent 


through all intervals, with the results of the 17 m to 11 m interval shown in Figs. 4.21 and 4.22. 


As expected, the estimated resonant frequencies are lower for Event 7, and damping is 
higher. The ability for Event 8 filters to so well model Event 7 events causes additional difficulty 
is interpreting these results - what are the “correct” values? Models created for Event 7, even using 
recursive techniques, sometimes do not map the data as “well” as the filter calculated from Event 
8, although results using the representative filter give physically satisfying results - we expect the 
damping to be greater for Event 7 since there was quite high peak acceleration and Arias intensity. 
Our interpretation is that the change in soil behavior between Event 7 and 8 is small enough that 
the low-strain filter captures the essence of the soil layer behavior. The statistical confidence of 
this estimate is illustrated in the graphs shown in fig. 4.23. This figure shows the 99% confidence 
interval in both the time and frequency domain. The slight change in soil properties between 


events is reflected in the Event 7 filter. 
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Figure 4.18 Comparison of the actual Event 7 surface layer output to the output estimated from 
the appropriate Event 8 filter 
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Figure 4.19 Comparison of the Event 7 E-W output to the Event 7 output modeled through the 
Event 8 N-S filter 
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Figure 4.20 Results of modeling the Event 8 N-S output through the Event 7 N-S filter, 6 m to the 
surface. 
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Figure 4.21 Results of modeling the Event 7 N-S output through the Event 8 N-S filter, 17 m to 


11 m. 
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Figure 4.22 Results of modeling the Event 8 N-S output through the Event 7 N-S filter, 17 m to 
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Figure 4.23 99% confidence interval for estimated output, Event 8, 6 m to the surface, in the time 
and frequency domain. 
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Figure 4.24 Comparison of pseudo-stiffness for DHA and DHB for Events 7 and 8. 
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4.4.1 Event 7 and 8 Soil-Structure Interaction 


Effects of soil-structure interaction (SSI) on the stiffness and damping estimates are evident 
from entries in Tables 4.1 and 4.2. As for Event 4, all effects of SSI are gone by a depth of 17 m, 
and are hardly evident at 11 m. These effects are evident in Fig. 4.24, which is a comparison of 
pseudo-stiffness for DHA and DHB for Events 7 and 8. The pseudo-stiffness is the same for DHA 
for both events, implying that the structure mitigated much of the strain-induced effects. The free- 
field stiffness is initially much greater than for DHA, and then approaches the level of the soil 
affected by SSI. For the second interval, 11 m to 6 m, there is only slight effects of SSI evident for 
Event 7, as shown in Fig. 4.25. For the deep two intervals, Figs. 4.26 and 4.27 show identical 


behavior for both arrays and both events. 
4.5 Event 16 


Event 16 was the highest magnitude (M, = 7.0) event recorded at Lotung, and also had the 
greatest Arias Intensity (3995 m/s). Pore water pressure was recorded for this event, and is 
included in Appendix C. A summary of response in comparison to estimated strain 1s given in Fig. 
4.28 (Zeghal and Elgamel, 1993). Although the strain was large enough that the estimated 
damping ratios were greater than for low energy events and the resonant frequency was greater (see 
Tables 4.1 and 4.2), it would be of tremendous interest to actually track softening of the soil 
throughout the event. An obvious indication of softening would be phase shift between the deepest 
layer and the surface record (e.g. Glaser, 1996). Figure 4.29 is a comparison of N-S velocity for 
the 47 m and surface layers. It is seen that the greatest effect of 47 m of soil is an “amplification” 
of velocity at the surface, which is evident for the 47 m to 17 m interval as well, Fig. 4.30. The 
only indication of phase offset is between about 24 s and 27 s, and it is not a systematic offset such 


as would occur from serious softening of the intervening soil layers. 


Another method of looking for changes.in the system through time is to use a segmentation 
algorithm which uses parallel recursive algorithms to find parametric changes within stated limits 
(e.g. Ljung, 1993; Andersson, 1985). Figure 4.31 presents the fundamental frequency and 


damping estimates for the 6 m to surface interval, and shows a change in the system that can be 
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Figure 4.25 Comparison of pseudo-stiffness for DHA and DHB for Events 7 and 8, 11 mto6 m. 
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Figure 4.26 Comparison of pseudo-stiffness for DHA and DHB for Events 7 and 8, 17 mto 11 m. 
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Event 7 and Event 8, Free—Field vs. SSI, 47 m to 17 m Pseudo-Stiffness 
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Figure 4.27 Comparison of pseudo-stiffness for DHA and DHB for Events 7 and 8, 47 mto 17 m. 
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Earthquake: LSST 16 Station: PN2-1 (Depth 6.30 m) 
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Figure 4.28 Pore water pressure buildup and shear strain history, Event 16, DHA, 6 m depth (from 
Zaghel et al., 1995). 
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Figure 4.29 Comparison of N-S velocity for the 47 m and surface layers, DHA, Event 16. 
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Figure 4.30 Comparison of N-S velocity for the 47 m and 17 m layers, DHA, Event 16. 
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Event 16, 6 m to Surface 
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Figure 4.31 Fundamental frequency and damping estimates for 6 m to surface interval, Event 16. 
Note change that can be attributed to softening at about 22 seconds into the temblor. 
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Figure 4.32 Frequency and damping estimates for 11 m to 6 m interval, Event 16. The variance 


needed indicate a system change is lower than for Fig. 4.31. 
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attributed to softening at about 22 seconds into the temblor, as per Zeghal et al. (1995). Figure 4.32 
is a Similar presentation for the 11 m to 6 m interval, but with the required change in parameters to 
indicate a system change held to closer tolerance. The change in system parameters correlates very 


nicely with the increase in pore pressure shown in Appendix C and Fig. 4.28. 


However, the change in parameters throughout the event is tenuous and sensitive to 
analysis method. The “average” values, although they reflect softening, can easily “overpower” 
the subtle temporal changes. An example is shown in Fig. 4.32, which shows the excellent fit for 
the 6 m to surface interval using a time-invariant model, and the frequency response for this 
system. The goodness-of-fit in the frequency domain is indicated by the 99% confidence interval 
given by the dashed lines. The question remains as to which model - time invariant or non- 


stationary - better “captures” the nature of soil behavior. 
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Event 16, DHB, Actual vs. Modeled Output, 6 m to Surface 
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Figure 4.32 Event 16, 6 m to surface, in the time and frequency domain using a time-invariant 
model, and the frequency response for this system. The goodness-of-fit in the 
frequency domain indicated by the 99% confidence interval in dashed lines. 
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CHAPTER 5 - CONCLUSIONS 


This report describes the process and results of using system identification to calculate soil 
system parameters for eight earthquakes of varying intensities recorded at the Lotung Large Scale 
Seismic Test Site in northeastern Taiwan. For this study the modal frequencies and damping ratios 
were calculated for depth intervals of 0-6 m, 6-11 m, 11-17 m, 17-47 m, 0-11 m and 11-47 m for 
events 3, 4, 7, 8, 9, 10, 12 and 16 with local magnitudes ranging from 4.5 to 7.0. The modal 
frequencies and damping ratios calculated are examined for the effect of local energy intensity and 
soil-structure interaction. Modal frequencies are seen to decrease with increasing intensity once a 
certain threshold of acceleration/intensity is reached. This result is consistent with the data 
obtained by other authors using different techniques. For the 0-6 m interval the decrease in 
frequency with event energy is less pronounced under a model containment structure than in the 
free field. This soil-structure effect is increasingly diminished with depth and absent by the 17- 
47m interval. Calculated damping values demonstrate an expected increase with input seismic 
energy. For the 0-6 m and 6-11 m intervals the damping values are higher under the model 
structure than in the free field. This distinction is completely missing in the 17-47 m results. The 
transition to non-linear behavior, while less pronounced with increasing depth, consistently occurs 


above a peak acceleration of 0.05 g or Arias Intensity of 100 m/sec. 


System identification was used to calculate modal frequencies and damping values for 
various intervals for eight earthquakes of a range of ground intensities. System identification is a 
technique that identifies the important parameters of a mechanical system by modeling the effect 
of the system on an input signal. Comparing the expected output predicted by an SI model with 
the actual output determines which modeled mechanical system is an accurate representation of the 
system being studied. This method was used to characterize the behavior of the soil at the Wildlife 
site in the Imperial Valley for the Elmore Ranch and Superstition Hills earthquakes of 1987 
(Glaser, 1996). The same author has also written a paper describing the potential advantages of 
parametric modeling (system identification) over spectral ratio descriptions of a system transfer 
function (Glaser, 1995). Essentially, parametric models give a more reliable estimate of signal 


spectral content when short, and/or non-stationary signals are being considered. A specific 
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parametric model, the ARMA (autoregressive-moving-average) model, is ideally suited for this 
problem because it can be derived from the equations of motion for an N degree of freedom (N- 
DOF) oscillator. The vibrating soil layer is the system of interest, the natural frequency and 


damping ratio information are contained in the system parameters. 


The results clearly indicate some non-linear response over the intervais studied. Evidence 
of a decrease in specific interval fundamental frequency and an accompanying general trend of 
increased damping with higher seismic energy are clear. Comparison of the results of this study 
with previous work considered with the inherent superiority of parametric modeling for transient 
and/or non-stationary time series such as earthquakes indicate that system identification is a more 
robust method for identifying fundamental frequencies and damping values for layers of earth 


materials when borehole information is available. 
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APPENDIX A: PROCEDURES FOR DATA PROCESSING 
A.1 Data Processing 
A.1.1 Cataloging 


The EPRI-supplied data consists of 10 pc-format floppy disks of data. The data supplies 
ground motions, structural motions, dynamic earth pressure, and ancillary data for 18 earthquake 
events (Tang and Tang, 1992). The data set consists of 2,103 individual files. Pore water pressure 


values are available for events 12, 16, and 17 only. 


The first step of the data processing was to enter the relevant data files into MATLAB and 
group them into corresponding event files. Files eventO1 through event18 were created. A 
simplified naming convention was developed for each data record. The first letter of the name is 
an a, v, or d, for acceleration, velocity, and displacement, respectively. The next three letters are 
dha or dhb for downhole array a or b. The following number refers to the depth in meters - 0, 6, 
11, 17, and 47. The second number, following the underscore refers to the event number. The last 


letter - e, n, or u, refers to accelerometer orientation - east-west, north-south, or up-down. 


An example of the naming convention is adha47_18u. This record is the acceleration time 
history at a depth of 47 m at downhole array a, event 18. The record gives the vertical ground 


motion for this location and event. 
A.1.2 Filtering, Resampling, and Integration 


The data as received from EPRI is in the form of raw acceleration records, which we 
processed following standard U.S.G.S. method (Converse and Brady, 1992). The acceleration 
records are digitized at a rate of 200 samples per second (s/s), for a bandwidth of 100 Hz. This 
relatively high Nyquist frequency causes the event records to be very long, 8,000 data points for 
each 40 s trace. In addition, there is little useful information for our study above 10-15 Hz, and 
this region would be very noisy. Resampling greatly reduces the amount of data to be later 


analyzed, and eliminates time-domain aliasing of the band-passed signals. It was therefore 
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decided to low-pass filter and resample the data at a rate of 25 s/s, yielding a Nyquist frequency of 
12.5 Hz. Previous work indicated that the information of interest would be contained in this band 


(Glaser, 1993, 1995a,b). 


The data was resampled using the resample algorithm from the MATLAB Signal 
Processing Toolbox (Krauss et al., 1994). The data is first low-pass filtered using a Kaiser- 
windowed linear-phase FIR filter using ten terms on either side of the timestep in the calculation. 
The low-pass filter is applied both forward and reverse to eliminate phase shift, and the data 
resampled at 25 Hz. The same process was carried out for every strong motion record as well as 


for the pore water pressure records. 


The acceleration records were then processed and integrated to provide velocity records, 
and the velocity records similarly reprocessed and integrated to yield displacement time histories. 
The acceleration record is first high-pass filtered at 0.15 Hz using a 4th order bi-directional 
Butterworth filter to remove dc offset and low frequency drift. A best-fit straight line from the 
arrival of the strong motion is then subtracted from the data, followed by the sample mean. The 
pre-arrival data is then set to zero. The processed acceleration data is now integrated using the 
trapezoidal method (Converse and Brady, 1992). The same steps are carried out on this newly 


formed velocity time history to yield the strong motion displacement. 
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APPENDIX B - MATLAB PROCEDURE TO INTEGRATE ACCELERATION RECORDS TO 
VELOCITY AND DISPLACEMENT. 


function [a,v,d]=vd(f,Lcut,dT,n,pre,zname) 

% function [a,v,d]=vd(f,Lcut,dT,n,pre) 

%o 

% detrends and filters the input acceleration, 

% and integrates twice to give velocity and displacement. 
% The results are plotted so they can be reviewed. 

v/a 

% f is the input acceleration vector 

% Leut is the low-end cutoff frequency in Hz. 

%o AT is the time step 

% n is the order of the Butterworth filter; 

% pre is the pre-event segment length to be zeroed. 

Jo 

% since the acceleration is filtered twice, the 

% effective order of the filter is double the value of n. 

% 

zip=(1:pre); 

zip=zeros(size(zip)); 

% : 

[b,c]=butter(n, Lcut*dT*2.0, ’high’); 

a=filtfilt(b,c,f); 

a=dtrend(a(:,1),1,pre); 

a=detrend(a); 

a=a-a(pre+1); 

a(1:pre)=zip; 

% 

v=inttrap(a,dT); 

v=filtfilt(b,c,v); 

v=dtrend(v(1,:)’,1,pre); 

v=detrend(v); 

v=v-v(pre+1); 

v(1:pre)=zip; 

% 

d=inttrap(v,dT); 

[b,c]=butter(n, (Lcut*dT*2)/1.0, ’high’); 
d=filtfilt(b,c,d)’; d=dtrend(d(:,1),1,pre); 

d=detrend(d); 

d=d-d(pre+1); 

d(1:pre)=zip; 

Jo 

time=(length(f)*dT)/1.0; 

t=[dT: dT: time]; 

temp=name(1:(length(name)-6)); 

clf 
subplot(3,1,1), plot(t,a(1:length(t))), title([temp,’: acceleration’]), grid on 
subplot(3,1,2), plot(t,v(1:length(t))), title([temp,’: velocity’]), grid on 
subplot(3,1,3), plot(t,d(1:length(t))), title({temp,’: displacement’]), grid on, xlabel(’Time (seconds)’) 
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APPENDIX C: PORE WATER PRESURE TIME HISTORIES FOR EVENTS 
12, 16, 17. 
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APPENDIX D: REVIEW OF DAMPING VALUES MEASURED IN THE FIELD AND 
LABORATORY 


D.1 Introduction 


In order to put the damping estimates made in this report in context, it is important to review 
previous estimates of soil damping. Measurements have been made in the laboratory as well as the 
field. Advantages of laboratory tests include full control over loading paths and strain, including 
large-sirain measurements. Laboratory tests will all yield at best approximate results since no one 
can run a laboratory test on an undisturbed loose sand. The preliminary reports from back- 
calculating earthquake response imply that the laboratory degradation curve might be too high at 
intermediate strains. In addition, the results from two independent methods (Chang et al., 1990; 
Abdel-Ghaffar and Scott, 1979) show that the customary hyperbolic shape of the laboratory 
damping curve might be incorrect, and actually is S-shaped. Laboratory tests also examine a very 


small volume of soil, giving a point estimate compared to the site of interest. 


Field testing avoids the disturbance issues associated with laboratory testing. In principle, 
geophysical field testing methods measure damping of undisturbed expanses of soil. However, 
the system damping thus calculated is the small strain damping, or D,;,, and is only valid for the 
elastic region of the soil (Dobry et al., 1982). The limitation is due to the inability to reliably impart 
strains into the soil much greater than 1x10°°. Therefore, it has been impossible to measure 


threshold strain, iP and the soil degradation curve, D/D,,;,, in situ. 


Attempts to input enough energy into the ground to cause intermediate to large strains have 
not been very successful. The coma of energy needed would destroy, a bore-hole, and would be 
destructive on the surface as well. There is also the problem of the transducers being in the near- 
field ¢f they are close to a source large enough to cause large strains in an immediate area. One 
exception was a project undertaken for the Nuclear Regulatory Commission (Shannon-Wilson, 


1976) where intermediate-to-large strains were input in a large scale field experiment. 
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The optimum situation would be the ability to make measurements during different 
magnitudes of earthquake excitement. In this case shear strain in the layers of interest, and stiffness 
(velocity), would be continually monitored. This approach is optimal since it allows 
nondestructive (in as much as an earthquake is nondestructive) evaluation of soil properties in the 
actual situation of interest. Since earthquakes can not be made-to-order, the chances of this 
situation happening are virtually nonexistent. The instrumentation would also be extremely 
difficult. However, use of inverse theory allows the soil parameters of interest to be calculated 
from attainable data — the ground motion records of the motion going into the layers of interest, 


and above the layer itself. 


The following sections present a non-exhaustive compilation of damping estimates made 
in the laboratory, in the field using geophysical methods, and in the field using strong motion 
excitation. The “Source” column in the accompanying tables refers to the publications tabulated 


at the end of this appendix. 
D.2 Laboratory Estimates of Soil Damping: 


Table D.1 presents a summary of some laboratory estimates of soil damping available in 
the literature. The values given are “average” representative values commonly accepted by the 
geotechnical community. The low-strain (< 0.001%) damping estimates for soils, whether 
cohesive or non-cohesive, lie in a very close range from about 1% to 5%. Evenata large strain of 
0.1%, damping is estimated to be 10 - 12% or less. At large strains, sand displays significantly 
higher damping than clay, as would be expected from application of Mindlin’s theory (Dobry et 
al., 1982). Details of the relationship between strain and damping for sands measured in the 
laboratory are given in Fig. D.1 (Seed et. al., 1986), while Fig. D.2 presents a similar relationship 
for normally consolidated clay (Dobry and Vucetic, 1987). 


A comparison of the strain dependent damping for sands and clayss based on an earlier 
compilation (Seed and Idriss, 1970) is shown in Fig. D.3, where it can be seen that damping 
increases with strain faster for sand than clay. A similar relationship is shown in Fig. D.4, which 


is based on laboratory measurements done by Japanese researchers (Kokusho, 1987). The results 
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from the extensive laboratory testing done for the Lotung project is given in Fig. D.5, where it is 
seen that damping is a bit higher than the Seed and Idriss approximation. While there has not been 
as much work done on gravelly soils, Fig. D6 implies that these soils have a damping response 
much the same as sand (Seed et al., 1986). These Figures show that there actually is significant 
scatter of soil damping above the threshold strain. 


Table D.1: Laboratory estimates of Soil Damping 


Ca 
Ca 


sand 1-4 10°3 15 
2-8 10°72 
7 107! 
sexe A 
silty sand 1077 15 


10° 
10°! 


silty sand, Wildlife site 


silty sand, Lotung site 
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Figure D.1 Damping ratios for sands based on laboratory measurements (Seed et al., 1986). 
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Figure D.2 Damping ratios for normally consolidated clays based on laboratory meaurements 
(Dobry and Vucetic, 1987). 
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Figure D.5 Summary of laboratory damping estimates for the Lotung site (Anderson, 1993). 
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D.3 Field Geophysical Estimates of Soil Damping: 


Representative values culled from the literature for field geophysical estimates of soil 
damping are presented in Table D.2. Note that these estimates were made for displacements well 
below the threshold strain, so there is no evidence of damping increase. The value of these esti- 
mates is that they were made on “undisturbed” soil, and should therefore be a more accurate 
measure of soil physical properties. Details on some of field tests in Japan are given in Table D.3 
(Kokusho, 1987). The estimates are more scattered than from the laboratory, with damping 
estimates higher than in the lab. This relationship is presented graphically in Figure D.7, which 
compares laboratory and field damping from Japanese researchers (Kokusho, 1987). Clay 
damping is seen to measure higher in the lab than the field, while for non-cohesive soils field 
measurements are much higher when measured in the field. Which values are “correct” is not 


known. 
D.4 Field Estimates of Soil Damping Made From Strong Motion: 


A detailed study of damping at the Lotung site was done by Zeghal et al. (1995). Summary 
plots of results from the field results are given in Figure D.8. Detailed results for Events 7, 12, and 
16 are given in Figures D.9 through D.17. The other major studies using field excitation are Glaser 
(1996) at the wildlife site, and Qi et al. (1989) at Mexico City using the random decrement method. 
Table D.4 presents summaries of the values estimated by these studies. It can be seen that the 
damping values estimated from strong motion excitation of field sites are higher than laboratory 


estimates. 
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Figure D.6 Comparison of damping ratios for gravelly soils and sands (Seed et al., 1986). 
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Table D.2: Field Geophysical Estimates of Soil Damping 
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Table D.3 Damping ratios of various soils measured in the field by Japanese researchers 
(Kokusho, 1987). 
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Table D.4: Estimates of Soil Damping Made From Strong Motion 
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(Zeghal et al., 1995). 
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Figure D.9 Equivalent shear moduli and damping ratios at Event 7, 6 m depth, DHB (Zeghal et 
al., 1995). 
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Figure D.10 Equivalent shear moduli and damping ratios at Event 7, 11 m depth, DHB (Zeghal 
et al., 1995). 
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Figure D.11 Equivalent shear moduli and damping ratios at Event 7, 17 m depth, DHB (Zeghal 
et al., 1995). 
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Figure D.12 Equivalent shear moduli and damping ratios at Event 12, 6 m depth, DHB (Zeghal 
éet.al., 1995): 
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Figure D.13 Equivalent shear moduli and damping ratios at Event 12, 11 m depth, DHB (Zeghal 
et al., 1995). 
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Figure D.14 Equivalent shear moduli and damping ratios at Event 12, 17 m depth, DHB (Zeghal 
et al., 1995). 
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Figure D.15 Equivalent shear moduli and damping ratios at Event 16, 6 m depth, DHB (Zeghal 
et al., 1995). 
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Figure D.16 Equivalent shear moduli and damping ratios at Event 16, 11 m depth, DHB (Zeghal 
et al, 1995): 
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Figure D.17 Equivalent shear moduli and damping ratios at Event 16, 17 m depth, DHB (Zeghal 
ét al, 1995). 
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APPENDIX E: COMPARISON OF ACTUAL INTERVAL OUTPUTS TO MODELED 
INTERVAL OUTPUTS FOR ALL ANALYZED EVENTS 
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Figure E.1 Modeled versus actual output, Event3, E - W, DHA, 0 to 6 m depth interval. 
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Figure E.2 Modeled versus actual output, Event3, E - W, DHA, 6 to 11 m depth interval. 
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Actual (solid) VS. RARX Simulated (dashed) Out for[441] Fit: 0.078 
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Figure E.3 Modeled versus actual output, Event3, E- W, DHA, 11 to 17 m depth interval. 
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Figure E.4 Modeled versus actual output, Event3, E- W, DHA, 17 to 47 m depth interval. 
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Figure E.S Modeled versus actual output, Event3, E - W, DHB, 0 to 6 m depth interval. 
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Figure E.6 Modeled versus actual output, Event3, E - W, DHB, 6 to 11 m depth interval. 
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Figure E.7 Modeled versus actual output, Event3, E - W, DHB, 11 to 17 m depth interval. 
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Figure E.8 Modeled versus actual output, Event3, E - W, DHB, 17 to 47 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.6654 
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Figure E.9 Modeled versus actual output, Event4, E - W, DHA, 0 to 6 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.9853 
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Figure E.10 Modeled versus actual output, Event4, E - W, DHA, 6 to 11 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.6284 
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Figure E.11 Modeled versus actual output, Event4, E- W, DHA, 11 to 17 m depth interval. 
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Figure E.12 Modeled versus actual output, Event4, E - W, DHA, 17 to 47 m depth interval. 
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Figure E.13 Modeled versus actual output, Event4, E - W, DHB, O to 6 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.8071 
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Figure E.14 Modeled versus actual output, Event4, E - W, DHB, 6 to 11 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.569 
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Figure E.15 Modeled versus actual output, Event7, E - W, DHA, 0 to 6 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.6192 
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Figure E.16 Modeled versus actual output, Event7, E- W, DHA, 6 to 11 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.3628 
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Figure E.17 Modeled versus actual output, Event7, E- W, DHA, 11 to 17 m depth interval. 
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Figure E.18 Modeled versus actual output, Event7, E - W, DHA, 17 to 47 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.6394 
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Figure E.19 Modeled versus actual output, Event7, E - W, DHB, 0 to 6 m depth interval. 
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Figure E.20 Modeled versus actual output, Event7, E - W, DHB, 6 to 11 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.3701 
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Figure E.21 Modeled versus actual output, Event7, E - W, DHB, 11 to 17 m depth interval. 
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Figure E.22 Modeled versus actual output, Event7, E - W, DHB, 17 to 47 m depth interval. 
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Figure E.23 Modeled versus actual output, Event8, E - W, DHA, 0 to 6 m depth interval. 
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Figure E.24 Modeled versus actual output, Event8, E - W, DHA, 6 to 11 m depth interval. 
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Figure E.25 Modeled versus actual output, Event3, E - W, DHA, 17 to 11 m depth interval. 
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Figure E.26 Modeled versus actual output, Event8, E - W, DHA, 47 to 17 m depth interval. 


160 


18 


Solid: Model output, Dashed: Measured output Fit: 0.1088 


Velocity (cm/sec) 


0 50 100 150 200 250 300 350 400 
Time (sec) 


Figure E.27 Modeled versus actual output, Event8, E - W, DHB, 6 to O m depth interval. 
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Actual (solid) VS. RARX Simulated (dashed) Out for [881] Fit: 0.2178 
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Figure E.28 Modeled versus actual output, Event8, E - W, DHB, 11 to 6 m depth interval. 
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Figure E.29 Modeled versus actual output, Event8, E - W, DHB, 17 to 11 m depth interval. 
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Figure E.30 Modeled versus actual output, Event8, E - W, DHB, 47 to 17 m depth interval. 


164 


18 


Velocity (cm/sec) 


Actual (solid) VS. RARX Simulated (dashed) Out for[221] Fit: 0.07001 


Time (sec) 


Figure E.31 Modeled versus actual output, Event9, E - W, DHA, 11 to 0 m depth interval. 
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Figure E.32 Modeled versus actual output, Event9, E - W, DHA, 47 to 11 m depth interval. 
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Figure E.33 Modeled versus actual output, Event9, E - W, DHB, 6 toOm depth interval. 
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Figure E.34 Modeled versus actual output, Event9, E - W, DHB, 11 to 6 m depth interval. 
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Figure E.35 Modeled versus actual output, Event9, E - W, DHB, 17 to 11 m depth interval. 
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Figure E.36 Modeled versus actual output, Event9, E - W, DHB, 47 to 17 m depth interval. 
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Figure E.37 Modeled versus actual output, Event10, E- W, DHA, 11 to 0 m depth interval. 
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Figure E.38 Modeled versus actual output, Eventl10, E - W, DHA, 47 to 11 m depth interval. 
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Velocity (cm/sec) 


0 50 100 150 200 250 300 
Time (sec) 


Figure E.39 Modeled versus actual output, Event10, E - W, DHB, 6 to 0 m depth interval. 
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Figure E.40 Modeled versus actual output, Event10, E - W, DHB, 11 to 6 m depth interval. 
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Figure E.41 Modeled versus actual output, Event10, E - W, DHB, 17 to 11 m depth interval. 
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Figure E.42 Modeled versus actual output, Eventl0, E- W, DHB, 47 to 17m depth interval. 


176 
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Figure E.43 Modeled versus actual output, Event12, E- W, DHA, 6 to O m depth interval. 
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Figure E.44 Modeled versus actual output, Eventl2, E- W, DHA, 11 to6m depth interval. 
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Figure E.45 Modeled versus actual output, Event12, E - W, DHA, 17 to 11 m depth interval. 
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Figure E.46 Modeled versus actual output, Eventl2, E - W, DHA, 47 to 17 m depth interval. 
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Figure E.47 Modeled versus actual output, Event12, E - W, DHB, 6 to 0 m depth interval. 
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Figure E.48 Modeled versus actual output, Event12, E- W, DHB, 11 to 6 m depth interval. 
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Figure E.49 Modeled versus actual output, Event12, E - W, DHB, 17 to 11 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.8077 
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Figure E.50 Modeled versus actual output, Eventl6, E - W, DHA, 6 to 0 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.5604 
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Figure E.51 Modeled versus actual output, Eventl6, E- W, DHA, 11 to 6 m depth interval. 
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Figure E.52 Modeled versus actual output, Event16, E - W, DHA, 17 to 11 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.9809 
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Figure E.53 Modeled versus actual output, Event16, E - W, DHA, 47 to 17 m depth interval. 
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Figure E.54 Modeled versus actual output, Eventl16, E - W, DHB, 6 to 0 m depth interval. 
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Solid: Model output, Dashed: Measured output Fit: 0.9194 
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Figure E.55 Modeled versus actual output, Event16, E - W, DHB, 11 to 6 m depth interval. 
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Figure E.56 Modeled versus actual output, Event16, E - W, DHB, 17 to 11 m depth interval. 
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